Development of a tight-binding potential for bcc-Zr. Application to the study of 
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We present a tight-binding potential based on the moment expansion of the density of states, 
which includes up to the fifth moment. The potential is fitted to bcc and hep Zr and it is applied 
to the computation of vibrational properties of bcc-Zr. In particular, we compute the isothermal 
elastic constants in the temperature range 1200-/^ <T < 2000K by means of standard Monte Carlo 
simulation techniques. The agreement with experimental results is satisfactory, especially in the case 
of the stability of the lattice with respect to the shear associated with C' . Ifowever, the temperature 
decrease of the Cauchy pressure is not reproduced. The T = OK phonon frequencies of bcc-Zr are 
also computed. The potential predicts several instabilities of the bcc structure, and a crossing of 
the longitudinal and transverse modes in the (001) direction. This is in agreement with recent ab 
^ ' initio calculations in Sc, Ti, Hf, and La. 

C/3 ■ 



PACS numbers; 62.20.Dc, 61.50.Lt 



I. INTRODUCTION 



^ ' A fundamental problem in Condensed Matter Physics is how to project microscopic interactions in many-body 
O ] systems onto a description in terms of a reasonably small number of degrees of freedom 0|. This is particularly 
important in computer simulation studies. Atomic-level simulations are, intrinsically, large-scale calculations and in 
spite of the huge improvement nowadays available in computing capacity, an efficient method to rapidly evaluate 
. energies which also treat forces in a physically realistic way, still remains a major difficulty to be solved. 
' The computation of solid properties always requires a particular choice for the interatomic potential. In the choice 
there is a compromise between physics and efficiency. Computational efhciency makes empirical or semiempirical 
potentials desirable, but at the same time one should require that the underlying physics behind the model potential 
is able to reproduce the properties of the system or at least those of interest. 
^-H • Several years ago Friedel |^ suggested that the starting point in the description of transition metals (TM) is a 
' band-picture with a strong d-character. In this sense, the remarkable parabola-like behaviour of the cohesive energy 
. and the bulk modulus exhibited by most TM as a function of the number of d-electrons ||^,^ clearly indicates that 
cohesion is mainly dominated by the d-states. This has motivated the development of many-body potentials based 
on the tight-binding (TB) approximation (For a review see Ref. Q). They are semiempirical in nature, which 
S makes them very appealing for computer simulation studies and at the same time they incorporate the band character 
^ of the metallic cohesion so that the attractive part of the interatomic energy turns out to be many body. 
(— I In the present work, we develop a TB potential based on the moment expansion of the density of states, which 

Q includes up to the fifth moment. It has been suggested that this is the lowest order needed to reproduce the general 
O ' trends in the elastic constants of hep and fee TM as a function of band filling . 

^ ' Zirconium, as is the case of other metals and alloys, is close packed at low temperatures, but undergoes a structural 
• I— I : phase transition of the Martensitic type to a bcc structure at higher temperatures ||,^,|lO| . Some features of this phase 
' transformation in Zr were previously studied by using the second moment approximation TB model pd] , p^ . Here, we 
^ , shall not focus on the structural phase transition itself but rather concentrate our interest on the clasti c properties of 



the high-temperature bcc-phase, which is accepted to be mainly stabilized by entropy effects [|13|,^ 15 1. In particular 
we calculate the temperature behaviour of the relevant elastic constants by using standard Monte Carlo simulation 
techniques. The results thus obtained are compared with the available experimental data. The agreement is rather 
satisfactory, but it does not allow us to draw conclusions about the reliability of the potential. This is provided by the 
computation of the phonon dispersion curves for the bcc-Zr at T = OK. Indeed, the comparison with recent ab initio 
calculations p^ ] is indicative that most of the fundamental physics governing the vibrational properties of bcc-Zr is 
contained in the interatomic potential model. 
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The paper is organized as follows. The next section is devoted to the construction of the interatomic potential. 
In Sec. Ill we describe the fitting procedure to Zr. In Sec. IV we present and discuss the results and in Sec. V 
conclusions are drawn. 



II. DEVELOPMENT OF THE INTERATOMIC POTENTIAL 



The interatomic potential is developed following the tight binding bond model (TBBM) by A. P. Sutton et al. 
l^jQ in the two-centre orthogonal approximation. The basis set of the TB Hamiltonian only includes d atomic 
orbitals and the crystal field interactions are neglected. 

The cohesive energy is decomposed into two terms, 

The term Epair is an empirical pair potential which stands for electrostatic and exchange-correlation interactions [ pT| , 
and the bond energy is 

Ebond = J2 n,{E){E-e,)dE, (2.2) 
■i 

where Ep is the Fermi energy, are the on-site Hamiltonian matrix elements in the atomic orbital representation, 
and ni{E) are the local densities of states (LDOS). 

In this TB model, the s-d hybridization, which is known to make a considerable contribution to the cohesive energy 
of TM |T^, is not treated appropriately. Nevertheless, it can be assumed that such contribution is included implicitly 
in the pair term, as in the paper by Girshick et al. pQ] , or either that this contribution should be proportional to the 
d band width and therefore it is included in the bond term. 

The condition of local charge neutrality is fulfilled by defining local Fermi energies. This is done through the 
relation, 

n,{E)dE ^ Nd. (2.3) 

where Nd is the number of d electrons per site and the atomic orbital energy is chosen to be the energy zero, = 
This method is equivalent to shifting the LDOS rigidly, whereas adjusting the on-site energies selfconsistently, as 
proposed by Sutton et al., means that the shifts of the LDOS are accompanied by a distortion of their shape. 

(n) 

The LDOS are constructed from their moments, /i^ , following the formalism of the recursion method of Haydock 
pl| as if the moments corresponded to an s band. In this way the computed LDOS are rotationally invariant [ p2| . 
That is, from the second to the fifth moments we compute the coefficients bi, ai, 62 and 02 of the recursion method. 
Since the coefficients a„, 6„ are convergent oscillating series their limit aoo, ^oo is estimated as Ooo = (oi + a2)/2 
and boo = [bi + &2)/2, and the coefficients a„, 6„ with n > 2 are assumed to be equal to a^o, boo which gives rise 
to the well known square root terminator of the continued fraction expansion of the diagonal elements of the Green 
function [ p4| . The integration of the LDOS in order to obtain the bond energy is performed numerically. 

To completely specify the bond energy we need to define the functional form of the off-diagonal Hamiltonian matrix 
elements in the atomic orbital representation, known as coupling strengths or hopping integrals. In the two-centre 
approximation the coupling strengths between sites i and j can be written as 

{ta\H\jl3) = {Kpir^^j/n.Wdda + Kpir.j /n,)Vdd. + <^Ur:j /n,)Vdds) R{n,) (2.4) 

where '^api'i^ij /^ij) is the angular dependence given by the symmetry of the d orbitals, as shown by Slater and Koster 
p5[ , and R{rij) is the dependence on distance of the dda, ddir and ddd bonds, which is assumed to be equal for all of 
them. The relative strength of the couplings Vdda ■ Vcid-n ■ VddS are the canonical values -6 : 4 : -1. 
For the radial dependence we take 

exp(— gr), r < ri 

R{r) — exp(— gr)(a + br + cr^){r2 — r), ri < r < r2 (2-5) 

r > r2. 

where a, b, and c are constants fixed by the condition that R{r) and its two first derivatives must be continuous at 
ri, which gives a = (3r^ — 3rir2 -I- ?'2)/(^2 — ri)^, b = (— 3ri -I- r2)/(r2 — ri)'^ and c = 1/(^2 — ri)^. By construction, 
R{r) is also continuous at the cutoff distance r2. 
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Finally the pair term of the cohesive energy is taken to be 



'pair ^ 




(2.6) 



with 



A{exp{-pir) 
A(exp(-pir) 



+ ^exp(-p2'')), ^ r < r 1 

+ (,cxp{—p2r)){d + br + cr^){r2 - r)"^ , ri < r < 7=2 

r > r2- 



(2.7) 







where a, 6, and c are constants fixed by the condition that Vpair{f) and its first two derivatives must be continuous at 
fi, which gives a = (6f^ — 4fif2 + r|)/(f2 — fi)*, b = 2(— 4fi + r2)/(?^2 — ^1)* and 5 = 3/ (f2 — ri)"'. By construction, 
Vpair (r) and its first derivative is also continuous at r2 ■ 

The model has, in principle, six fitting parameters: A, ^, pi, p2, VddS and g. Nevertheless, the number of electrons, 
Nd, and the cutoff distances ri, r2, ri, and r2 are also used as fitting parameters although they are not completely 
free since they are constrained by their physical character. 



In this section we discuss the procedure used to fit the model parameters to the zero temperature properties of Zr, 
both in the hep and the unstable bcc phases. 

Since the bcc structure is only stable at T > I135K, the values of some properties used in the fitting procedure are 
extrapolated from the high temperature experimental data, whereas others are obtained from ab initio calculations 
at T = OK. Most of the ab initio results have been calculated in this work using the WIEN97 code [Q, which follows 
the Full Potential Linearized Augmented Plane Wave method within density functional theory. All the computed 
quantities are obtained in the generalized gradient approximation (GGA) by Perdew, Burke and Ernzerhof P7| for 
the exchange-correlation potential. 

The properties of bcc-Zr that we have considered in the fitting procedure are: the cohesive energy, the lattice 
parameter, the elastic constants and the unrelaxed vacancy formation energy. Moreover, we have also taken some 
properties of hcp-Zr at T = OK into account, mainly its elastic constants, lattice parameter, the c/a ratio, and the 
energy difference between both structures. 



The coupling strengths are essentially short-range functions. Therefore, we impose that they fall to zero between 
the second and third nearest neighbours and choose r2 — 1.24abcc- For the matching point of the tail of the function 
we chose ri = 0.92ai,cc, which lies between the first and second nearest neighbours. Notice that in this way the tail of 
the coupling strengths interferes with the fitting procedure. The values chosen for these parameters are those which 
allow a better reproduction of the properties of hcp-Zr. 

An important property of the elastic constants exploited in the fitting procedure is that the Cauchy pressure given by 
an interatomic potential which satisfies the mechanical equilibrium conditions only depends on the many-body term of 
the potential, that is, the Cauchy pressure is independent of the pair term |20[ |. Considering both the bcc and the hep 
structures we have three Cauchy pressures to deal with. One of these, the Cauchy pressure C12 — Cqq of the hep lattice 
is affected by an internal relaxation of the lattice under strain. In order to avoid the determination of the internal 
relaxation and its effect on the Cauchy pressure during the fitting procedure, as the test value for this quantity, we 
use only its homogeneous contribution (without internal relaxation) obtained from the experimental Cauchy pressure 
by subtracting the inhomogeneous contribution determined from ab initio calculations. In the hep lattice the elastic 
constants affected by internal relaxations arc Cn, C12 and Cee — (Cn — Ci2)/2, and the inhomogeneous contribution 
to C12 is the same as that of Cn but with the opposite sign. Therefore, we need only to estimate the inhomogeneous 
contribution to Cn which is / ~ 7 GPa Using the experimental values of the elastic constants together with the 
ab initio results, the three Cauchy pressures are, 



III. FITTING TO ZR 



A. Fitting of the bond energy 




(3.1) 
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Another property which strongly depends on the many-body term is the unrelaxed vacancy formation energy, 
{Ey)^'^'^°^^ , since the contribution of the pair potential term to this quantity is equal to its contribution to the 
cohesive energy. For the present interatomic potential, in the bcc lattice we get the following approximate relation, 

{Ey)^^^°^^ — —Epair — 0.52Ebond- (3.2) 



Taking into account (2.1) we therefore get, 

Ehond ^ {{Elr^''^'^ + E,,h)/0A8 ~ -8.0 eV. (3.3) 

This result is consistent with the value predicted by the renormalized-atom method of Gelatt et al. (lo) provided 
that the contribution to the cohesive energy of the s-d hybridisation is assumed to be contained in the bond term 

(-^d band broadening "t" -£'s—d hybridization — 7.8eV). 

From the literature we obtain physical values for the number of electrons, Nd — 2.536 (29j, band width, W=7.8 eV 



and the derivative of the band width with respect to the atomic volume, —3ild\nW/dfl = 3.97 30,3l|], which is 
related to the parameter q. In order to extract the value of the parameter q from this relation we approximate the 
band width as (x that is, proportional to the square root of the second moment of the density of states 

and using the radial form of the coupling strengths defined in Eq. (|]^), we obtain q ~ 4.36a^^^. In the present model 
the band width is W = Etop — Ebottom = 46oo which is proportional to the parameter VddS- Therefore we can use W 
instead of VddS as the fitting parameter without loss of generality. 

Using these physical values for the parametres Nd, W and g, in the bcc lattice the interatomic potential gives a value 
of E iiond — —7.5 eV, which is consistent with the value predicted by the unrelaxed vacancy formation energy (Eq. 
( |3.3|) ). Nevert heles s, the Cauchy pressure obtained with these parameters, (C12 — Ci4)bcc — 104 GPa, is unacceptably 
high (see Eq. (3.1)). We have, therefore, to find an alternative way of fitting them. We shall proceed by paying more 



attention to the directly observable physical quantities such as cohesive energy, lattice parameter, elastic constants 
and vacancy formation energy rather than to the physical interpretation of the parameters. 

For a given value of Nd, the parameters W and q are determined numerically to reproduce the Cauchy pressure of 
the bcc lattice given in Eq. (|]|) and the bond energy, E^ond = — S.OeV (Eq. ( |3.3[ )). We repeat this fitting procedure 
for different values of Nd and compute the Cauchy pressures of the hep lattice. The results are given in Table |. 
These Cauchy pressures are computed using the lattice parameter and c/a ratio of the hep lattice at which this 
structure is stabilized by the interatomic potential. For this purpose we use an arbitrary pair contribution fitted to 
the lattice parameter, cohesive energy and elastic constants of the bcc structure. In all cases, aucp and c/a are close 
to the experimental and ideal values respectively, although we observe a tendency of c/a to increase as the number of 
electrons is increased. The value of Nd which allows reproduction of the three Cauchy pressures simultaneously (Eq. 
( |3.1[ )) is Nd = 1.45, which is substantially lower than its physical value, Nd = 2.536 |^^. Nevertheless, there are two 
additional reasons for adopting the value of Nd — 1.45, i) TB studies on the relative stability of hep and fee lattices 
as a function of Nd have shown that the hep lattice is only stable for Nd < 2 (in the range 0.5 < Nd < 4.5), and 
ii) the low value of the c/a ratio in group IV hep metals can be justified in terms of the different contribution to the 
bond energy from nearest neighbours located on the basal plane and those located off it only if Nd < 2 jl^,^ ■ 

The radial dependence of the coupling strengths is shown in Fig. |l|. 



B. Fitting of the pair potential 

The cutoff distance of the pair potential is chosen to be r2 = 2.211abcc which lies between the 7th and 8th nearest 
neighbours in the bcc structure. This value has been adjusted in order to reproduce the properties of the hep lattice. 
The value of fi is almost irrelevant since it does not really affect the functional form of the pair potential. We chose 
fi — l.2abcc- For given values of the parameters pi, p2 and ^, the parameter A is determined from the mechanical 
equilibrium condition of the bcc lattice 



dEcoh 



da 



0. (3.4) 



The parameters pi, p2 and ^ are given by the cohesive energy, and elastic constants Cn and C12 of the bcc lattice. 
Notice that the elastic constant C44 is automatically reproduced by fitting C12 since the Cauchy pressure is fixed by 
the bond energy. 

The functional form of the pair potential is shown in Fig. |^. 
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The values of the fitted parameters are given in Table || and the properties of bee and hep Zr used for the fit are 
given in Tables Ul and IV respectively. In the fitting procedure we have priorized the reproduction of the properties 
of bcc-Zr since the interatomic potential will be used to study the vibrational properties of this structure. This have 
led to different accuracy in the fitted quantities of bcc and hep Zr (see Tables III and |^). 



IV. CALCULATION OF THE VIBRATIONAL PROPERTIES OF BCC-ZR 



In this section we use the interatomic potential already obtained to compute the phonon dispersion curves of the 
bcc lattice at T = OK, and the elastic constants in the temperature range 1200K < T < 2000K, where experimentally, 
the bcc structure is stable. 



A. Phonon frequencies of bcc-Zr at zero temperature 



In order to obtain the phonon frequencies of bcc-Zr, we first compute the interatomic force constants by means of 
standard numerical differentiation. It is worth pointing out that the interatomic potential gives rise to long-range 
force constants. In particular, the range of the pair potential contribution is the range of the potential, which in 
the present case is up to the 7th nearest neighbours. The range of the bond term contribution is much larger. For 
coupling strengths extending up to the second nearest neighbours and including up to the fifth moment of the LDOS 
in the computation of the bond energy, the range of the force constants extends up to the 22th nearest neighbours. 
This is due to the many-body character of the bond energy together with the dependence of the high order moments 
on the position of distant atoms. 

In Table ^ we show the results for the computed force constants. From these values we construct the dynamical 
matrix and compute the phonon frequencies in the harmonic approximation along the high symmetry directions of 
the reciprocal space (Fig. ||). Since the interatomic potential is fitted to the T = OK elastic constants of bcc-Zr, the 
slope of the phonon dispersion curves around the T point is expected to be correct. There are several features of the 
phonon dispersion curves that deserve special comment: 

1. The whole Tl(^ ^ 0) branch is unstable. 

2. The Tl(^ ^ 2^) branch has a positive slope around the F point (consistent with the associated combination of 
elastic constants), but it rapidly softens and becomes unstable. Before matching the Tl(^ ^ 0) branch at the N 
point, it becomes stable again at about ^ = 1/3. 

3. The softening of the L(^ ^ ^) branch around ^ = 2/3 observed experimentally at high temperature is, at 
zero temperature, an instability which gives rise to the w phase. The frequency of the ^ = 2/3 phonon is 
approximately zero, as was obtained by K.M. Ho et al. from ab initio calculations, and the minimum of the 
branch is at about — 17/24 

4. In the (001) direction there is a crossing between the longitudinal and transverse branches. 

Several of these features, not observed experimentally at high temperature, were obtained from ab initio calculations 
in Sc, Ti, Hf, and La at T = OK In these materials the whole Tl(^ ^ 0) branch has imaginary frequencies (the 
phonons in the (112) direction have not been computed). The instability around the ^ = 2/3 L(^ ^ ^) mode is also 
predicted, although the minimum of the branch is not located at ^ = 17/24 but at f — 7/12 for all elements (except 
Sc, which has the minimum at ^ = 2/3). Finally, the crossing of the longitudinal and transverse branches in the (100) 
direction is also observed. 

Moreover, in these materials the T(^ ^ ^) branch has imaginary frequencies around the F point. The slope of this 
branch is given by the elastic constant Cn — C12 -I- C44, which in Zr is positive, and thus this feature cannot be 
observed. 

From the force constants we also compute the elastic constants following the method of long waves |Q , and recover 
the values obtained by means of homogeneous deformations. This can be used as a test of the internal consistency of 
the interatomic potential |l^. 
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B. High temperature elastic constants of bcc-Zr 



The bcc-Zr high-temperature elastic constants are obtained from Monte Carlo (MC) simulations in the canonical 
ensemble, {T,V,N), using the atomic volume obtained from MC simulations in the isobaric-isothermal ensemble, 
{T,P,N), at zero pressure. The theoretical background of such simulations can be found in the work by I. R. 
McDonald ||. 

The second-order isothermal elastic constants are computed using the fluctuation formula 
„T I d^F 1 / d^Ec 



V , V , V ^ (4.1) 

1 (/ dEcoh dEcoh\ dEcoh\ dEoohW NkeT , ^ ' 



ksTV I \ dei-j deu / \ deij / \ deu \ V 



where F is the Helmholtz free energy, are the elastic strains, V is the total volume, is the Boltzmann constant, 
T is the temperature, and N is the number of atoms. 

The derivatives of the cohesive energy with respect to the elastic strains are computed numerically. Since the second 
derivative of both the pair potential and the coupling strengths is not continuous at the cutoff distance, the numerical 
derivative must be calculated using the same radial and pair functions in both the strained and the unstrained state 
of the lattice, for each of the atomic pairs. That is, if for a given pair of atoms the distance, r, is ri < r < r2, the 
radial function used for the computation of the bond energy of both the strained and unstrained lattice will be that 
defined in this region, regardless of the fact that in the strained lattice we may have r > r2. 

The simulations are performed ona4x4x4x2 = 128-site bcc lattice with periodic boundary conditions. The 
attempted configurational changes are single atom displacements and after each N of these attempts {— 1 Monte 
Carlo step (MCS)), in the isobaric-isothermal ensemble a volume change is also proposed. The simulations are 10^ 
and 5 • 10* MCS long in the isobaric-isothermal and canonical ensemble respectively. 

Due to the many-body character of the bond energy, the change in the total cohesive energy due to a single- 
atom movement involves the recalculation of the contribution to bond energy of about 65-110 atoms, depending on 
temperature. Nevertheless, this computation can be highly optimized and only requires about six times the CPU 
time needed to compute the contribution to the bond energy of a single atom. 

In Fig. H we show the computed lattice parameter vs. temperature. Although at temperatures above 12GGK the 
computed lattice parameter is about 0.6% smaller than the experimental value, the thermal expansion coefficient 
(slope) is reproduced to great accuracy (3 = 1/V (dV/dT) = 3.0 • IQ-'^K-^ {ficxp = 2.8 • IQ-^K-^ |^). Notice that 
the linear extrapolation of the computed lattice parameter to zero temperature does not match the fitted value. This 
is because at low temperature the thermal expansion given by the interatomic potential is strongly nonlinear. This is 
probably related to the fact that the bcc structure is unstable, although this behavior is not observed when comparing 



the high temperature experimental results | |37[ | with the zero temperature ab initio calculations (see Table III ) . 

In Fig. 1^ we show the temperature dependence of the relevant elastic constants obtained from Monte Carlo 
simulations. The main success of the present interatomic potential is that it renders a positive C at high temperature. 
Moreover, the value predicted for the whole set of elastic constants, Cn, C12, and C44 is rather accurate. The main 
failure is that it is unable to reproduce the large value of C" observed experimentally. This failure comes mainly 
from the values of the elastic constant C12 which experiments have shown to decrease strongly with temperature. 
This marked decrease of C12, together with the nearly constant behaviour of C44, means that the Cauchy pressure 
decreases with temperature. We were unable to reproduce such behaviour. In several tests during the parametrization 
of the interatomic potential we always found a Cauchy pressure nearly independent on temperature. 



In Table VI we show separately the different contributions to the elastic constants: the Born term, the fluctuation 
term, and the kinetic term. We get the rather unusual result that the fluctuation term of C12 is nearly zero. This 
means that the strains in different directions are uncorrelated in the simulations 

dEcoh dEcoh \ ^ , / '^^coh \ / dEcoh \ .^2-) 
de^x deyy / ~ \ de^^ / \ deyy / ' 

Since the contribution of the fluctuation term to C12 is usually negative and the interatomic potential is unable 
to reproduce the low value of C12 observed experimentally, we conclude that this lack of correlation given by the 
interatomic potential is possibly unphysical. 
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V. DISCUSSION AND CONCLUSIONS 



In the present paper we have developed a TB interatomic potential suitable for the study of the vibrational properties 
of bee Zr. The interatomic potential has been fitted to the T = OK properties of Zr in both the hep and the bcc 
structures. Although among the vibrational properties, only the elastic constants are used in the fitting procedure, 
the TB potential shows a remarkable capacity of predicting the T = OK phonon frequencies of the bcc structure 
along the high symmetry directions studied. As regards the high-temperature elastic constants, the general trends are 
reproduced, especially the stability of the bcc structure with respect to the shear associated with the elastic constant 
C". Nevertheless, the interatomic potential is unable to reproduce the temperature decrease of the Cauchy pressure. 

The reliability of the experimental values of the high-temperature elastic constants should, however, be questioned. 
The elastic constants cannot be obtained from measurements of the velocity of acoustic waves in the material because 
the temperature at which the bcc phase is stable is too high. Helming et al. therefore derived the elastic constants 
from the force constants obtained from a fit to the phonon dispersion relations. In order to keep the number of force 
constants reasonably small, in the fitting procedure they impose that the range is up to the fifth neighbour shell, 
which is rather short. On the other hand, the elastic constants obtained depend critically on the frequencies of the 
phonons close to the F point of the BriouUin zone, and thus have large error bars. We have tried to derive the elastic 
constants from the phonon frequencies of Helming et al. |^ , but we only were able to reproduce the elastic constant 
C to any accuracy. The values of all the other clastic constants strongly depend on the phonons considered and the 
method of fitting. 

In spite of the remarkable success of the interatomic potential in reproducing the T = OK phonon frequencies, 
we should mention the difhculties we have found during the fitting procedure, and discuss which features of the TB 
potential are expected to correctly reproduce the physics of the material and which are not. 

The first point concerns the range of the hopping integrals, which in fact is too small. In the hep lattice they should 
fall between the second and third nearest neighbours, at least, and only the nearest neighbours are taken into account. 
This leads to a bond energy contribution in the hep lattice that is smaller than in the bcc lattice, which is just the 
opposite to the expected result. This fact is compensated by the pair potential contribution to give the correct energy 
difference between both structures, but it is still clearly reflected in the unrelaxed vacancy formation energies. 

The reason for choosing such a low value for the cutoff distance of the hopping integrals is computational convenience. 
The range of the hopping integrals in the bcc lattice is up to second nearest neighbours (14 atoms involved). This 
means that in the perfect lattice at T = OK the number of neighbours involved in the computation of the moments is 
64 (up to the sixth coordination shell). Nevertheless, at T = 2000K the atoms are far from their equilibrium positions 
and this rapidly increases the number of neighbours involved in the computation of the moments up to ~ 110. In order 
to correctly compute the moments we must therefore take into account the neighbours up to the 13th coordination 
shell (258 atoms). An increase in the cutoff distance of the hopping integrals involves an increase of the number of 
neighbours to be taken into account in the computation of the moments, and thus, we decided to choose the lowest 
possible value which allowed us to obtain physically reasonable results. 

The second point concerns the pair potential contribution. During the fitting procedure we observed that the 
capacity of the interatomic potential to simultaneously reproduce the properties of both the hep and the bcc structure 
is strongly dependent on the range of the pair potential. In other words, if we take a different range to that used in 
this paper, the results are rather worse. This is indicative that the geometry of the different coordination shells has 
an important effect on the elastic constants. Moreover, although at the cutoff distance the pair potential and its first 
derivative are continuous, the decay to zero is still sharp, and the contribution to the elastic constants by the last 
coordination shell is unphysically high. 

Finally, we should mention that the s-d hybridization, which is not explicitly included in the TB potential, has an 
important contribution to the cohesive energy (about -2eV jl^). We have considered only d atomic orbitals in the 
basis set in order to minimize computation time and the complexity of the TB potential. 

The problems encountered when trying to use the physical values for the quantities Nd, W, and q have already 
been discussed. Nevertheless, the treatment of these quantities as fitting parameters gives enough flexibility to the 
interatomic potential to reproduce to reasonable accuracy all the magnitudes described in the paper. 

A significant improvement in the behaviour of the elastic constants requires a better determination of the Fermi 
energy, together with a more detailed description of the DOS, especially around this point. Nevertheless, the inclusion 
of high order moments into the interatomic potential is computationally very expensive. 
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TABLE 1. Cauchy pressures of hep Zr for different values of Nd- The parameter q and the band width W are also shown. 



Nd 


1 («bcc) 


W (eV) 


(C7i2 - Cm)tr'' (GPa) 


(Ci3 - C44)hcp (GPa) 


1.35 


3.730 


13.8 


1 


28 


1.40 


3.680 


13.4 


5 


30 


1.45 


3.634 


13.0 


9 


32 


1.50 


3.592 


12.6 


13 


34 


1.60 


3.498 


12.0 


21 


37 


1.80 


3.342 


10.9 


34 


44 


2.00 


3.201 


10.1 


44 


49 


2.20 


3.099 


9.2 


53 


53 


2.40 


3.012 


8.9 


57 


55 


2.60 


2.953 


8.5 


62 


57 



TABLE 11. Parameters of the interatomic potential and their physical values from ab initio calculations. 



Parameter 


Interatomic potential 


Ab initio 




1.45 


2.536 |2i 


1 («6cc) 


3.634 




W (eV) 


13.0 




ri {a^cc) 


0.92 




r2 (abcc) 


1.24 




A (eV) 


640.476022 






8.45462323 




^ (dimensionless) 


-1.80733296E - 05 




P2 (abcc) 


-0.961306397 




fi (abcc) 


1.20 




f2 (abcc) 


2.211 





TABLE III. Properties of bcc-Zr used for the fit. 





Interatomic potential 


Experiment / Ab initio 


a (nm) 
Ecoh (eV) 
Cii (GPa) 
Ci2 (GPa) 
C44 (GPa) 
C" (GPa) 


0.3574 
-6.15 
81.7 
93.4 
36 
-5.8 
2.32 


0.3574 iTf, ( 
~ -6.1 
81.7 1 
93.41 
36 1 
-5.8 
2.30 1 


L3580 ^ 
50 

2^1 
21 

sf 

|| 

3f 


TABLE IV. Properties of hcp-Zr used for the fit. 



Interatomic potential Experiment / Ab initio 

a (nm) 

c/a 

{Ecoh)bcc — {Eco (eV) 
^homog. ^gp^^ 

^homog. (.Qp^^ 

Ci3 (GPa) 
C44 (GPa) 
C33 (GPa) 
Cfis""''- (GPa) 



0.3196 

1.6284 

0.044 

161.8 

60.1 

68.2 

36.6 

179.5 

50.8 

2.44 



0.3229 M 
1.592 g 
0.04 HTi, 0.09 




TABLE V. Force constants of bcc-Zr obtained from the tight binding interatomic potential (in 10 ^ N/m). 



shell 


coord. X 2 


XX 


yy 


zz 


yz 


xz 


xy 


1 


111 


-5495.93 


-5495.93 


-5495.93 


-13410.85 


-13410.85 


-13410.85 


2 


200 


4320.75 


-6704.46 


-6704.46 











3 


220 


-454.85 


-454.85 


-620.49 








-932.34 


4 


311 


-2418.67 


-95.44 


-95.44 


-410.94 


-437.91 


-437.91 


5 


222 


592.90 


592.90 


592.90 


306.61 


306.61 


306.61 


6 


400 


584.12 


-42.25 


-42.25 











7 


331 


957.30 


957.30 


191.96 


346.44 


346.44 


934.45 


8 


420 


-73.56 


21.17 


14.13 








-55.33 


9 


422 


-122.99 


7.54 


7.54 


-7.77 


-70.19 


-70.19 


10 


333 


-98.35 


-98.35 


-98.35 


-107.07 


-107.07 


-107.07 


10 


511 


-82.20 


10.22 


10.22 


-3.17 


-14.81 


-14.81 


11 


440 


-34.35 


-34.35 


19.05 








-44.28 


12 


531 


-29.68 


-10.22 


3.53 


-3.35 


-8.72 


-23.40 


13 


442 


-14.60 


-14.60 


-3.05 


-8.83 


-8.83 


-17.23 


13 


600 


-64.36 


6.10 


6.10 











14 


620 


-27.97 


-3.05 


3.05 








-12.88 


15 


533 


-8.94 


-3.53 


-3.53 


-3.53 


-6.23 


-6.23 


16 


622 


-13.71 


-1.52 


-1.52 


-1.52 


-6.30 


-6.30 


17 


444 


-6.47 


-6.47 


-6.47 


-6.47 


-6.47 


-6.47 


18 


551 




















18 


711 


-12.22 











-2.66 


-2.66 


19 


640 




















20 


642 




















21 


553 




















21 


731 




















22 


800 


-4.67 


















TABLE VI. Born, fluctuation, and kinetic contributions to the second-order isothermal elastic constants (in GPa) obtained 
from the Monte Carlo simulations in the canonical ensemble at different temperatures (in K) and zero pressure. 



T 




'^12 


/^Born 
O44 


Qfluct. 


y-~lfluct. 








O44 


1188 


126.8 


92.6 


48.6 


-38.3 


3.4 


-16.3 


1.4 





0.7 


1300 


128.2 


92.4 


48.4 


-39.5 


3.7 


-16.4 


1.5 





0.8 


1400 


129.8 


92.2 


48.2 


-37.9 


1.7 


-16.5 


1.7 





0.8 


1483 


131.0 


91.8 


48.0 


-39.7 


1.9 


-16.8 


1.7 





0.9 


1600 


131.7 


91.7 


47.7 


-40.0 


0.9 


-17.0 


1.9 





0.9 


1700 


132.5 


91.5 


47.5 


-40.1 


-0.4 


-17.3 


2.0 





1.0 


1800 


133.1 


91.1 


47.2 


-40.4 


-1.4 


-17.2 


2.1 





1.1 


1883 


1.33.5 


90.7 


46.8 


-40.8 


-2.2 


-17.9 


2.2 





1.1 


2(KK) 


i:',r,X} 


<)().8 


Ki.l) 


-li.l 


-:!.<) 


-17.1 


2.:i 





1.2 



10 




11 




12 



0.365 



0.364 



a 



0.358 




1100 1300 1500 1700 1900 2100 

Temperature (K) 



FIG. 3. Lattice parameter of bcc-Zr at zero pressure vs. temperature. The solid line is the experimental result and the circles 
are the computed values. The statistical error is denoted by the size of the circles. 
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1100 1300 1500 1700 1900 2100 

Temperature (K) 



FIG. 4. Isothermal second-order elastic constants of bcc-Zr vs. temperature. Empty symbols are the computed values, and 
filled symbols are the experimental results. The horizontal dashed line emphasizes the change of the sign of the elastic constant 
C and the solid lines are guides to the eye. The elastic constants are Cn (circles), C12 (squares), C44 (triangles) and C 
(diamonds) . The statistical error of the computed values is smaller than the size of the symbols 
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